
clear;
n = 64;
N = n*n;

%set projection num
 L = 180;
 
%measure matrix
[matrix] = getSampMatrix(n,L);
A = double(matrix);

%static sample
% b = A * samplePIC(:);
% bback = b;

theta = (180/L ):(180/L):180;
%theta = thc*180/pi;
%[R] = radon(samplePIC,theta);
%Rback = R;
%xp = iradon(R,theta,n);
%xpStatic = xp(:);

%dyna sample
scale = n/64;
cycleNum = 5;
cycleLen = floor(L/cycleNum);
halfCycleLen = floor(cycleLen/2);
MaxDis = 4*scale;
speed = MaxDis/halfCycleLen;
ObjSize = 4*scale;
if halfCycleLen == 0
    disp('need more projection')
    return;
end

for i = 1:L
    j = mod(i,cycleLen);
    pos = j - halfCycleLen;
    dis = pos * speed;
    dis = floor(dis);
    a(j+1) = dis;
   
    center = n/2;
 
    samplePIC = SampPIC(center+dis,center+ObjSize+dis,center,center+ObjSize,n);

    Rtemp = radon(samplePIC,theta(i),64);
    R(:,i) = Rtemp;
    Ttheta(j+1,ceil(i/cycleLen)) = i;
    
     btemp = A * samplePIC(:);
     %detVal((i-1)*n+1:i*n,1) = dynaSamp(L,i,samplePIC(:),n);
     %detVal((i-1)*n+1:i*n,1) = dynaSamp(L,i,samplePIC(:),n);
   %-----------------add noise----
%    noiseLevel = max(btemp) * 0.1;
%     noise =   randn(size(btemp));
%     btemp = btemp +    noiseLevel * (noise/norm(noise));
	%-------------------------------
        proj(j+1,(ceil(i/cycleLen)-1)*n+1:(ceil(i/cycleLen))*n) = (i-1)*n+1:i*n;
        bcycle(j+1,(ceil(i/cycleLen)-1)*n+1:(ceil(i/cycleLen))*n) = btemp((i-1)*n+1:i*n);
     b((i-1)*n+1:i*n,:)= btemp((i-1)*n+1:i*n,:);       
end

xpDyna = iradon(R,theta,'spline','Shepp-Logan',n);



%Phase = [1:3,10,11,30,31];
Phase = [1:5];

%
detVal =   bcycle(Phase,:)';
detVal = detVal(:);
PrjSeq = proj(Phase,:)';
PrjSeq = PrjSeq(:);
thetaSeq = Ttheta(Phase,:)';
thetaSeq = thetaSeq(:);
%detVal  == A(PrjSeq,:) * xrec ;   

totalIter = 200;
artIter = 5;
GDIter = 20;
tol = 0.01;
%x0 = xpDyna(:);
x0 = zeros(length(xpDyna(:)),1);
error = 0;
%for xtrue
dis = -4;
xtrue = SampPIC(n/2+dis,n/2+ObjSize+dis,n/2,n/2+ObjSize,n);
xtrue = xtrue(:);

iterError = 0;
ObjV = 0;
lamda = 0;
%x0 = zeros(length(x0),1);
xlast = x0;

Aeq = A(PrjSeq,:);
beq = detVal;


tic;
for i=1:totalIter
    
    disp('totalIter');
    disp(i);
    
    
    %x0 = reconART(detVal,A(PrjSeq,:), tol, x0(:),artIter*length(detVal)); 
    
    x0 = reconART(detVal,Aeq, tol, x0(:),artIter*length(beq));
    %for plot
    iterError((i-1)*2+1) = norm(x0 - samplePIC(:));
    %iterError((i-1)*2+1) = norm(x0 - samplePIC(:))/norm(x0);
    xart = x0;
    %ObjV((i-1)*2+1) = ObjFun(x0,n,xart);

    
   % x0 = GradDes(xpDyna(:),lamda,x0(:),n,GDIter,tol,detVal,A(PrjSeq,:));
   x0 = GradDes(xpDyna(:),lamda,x0(:),n,GDIter,tol,x0(:));
    
    %--------------hard threshold-------------------
%         [temp_ht index_ht]= sort(x0,'descend');
%         ht_threshold = 838;
%         for ht=ht_threshold:length(x0)
%             x0(index_ht(ht))=0;
%         end
    
    %------------------------------------------------     
     iterError((i-1)*2+2) = norm(x0 - samplePIC(:));
     %iterError((i-1)*2+2) = norm(x0 - samplePIC(:))/norm(x0);
     %ObjV((i-1)*2+2) = ObjFun(x0,n,xart);
     xlast = x0;
    
%     if(iterError(i)<1)
%         break;
%     end
    
end

 timeDura = toc;



xresult = x0;

disp('finished');

resultPIC = figure;
subplot(3,1,1);

showPIC(xpDyna(:),n,n);
title('prior image');

subplot(3,1,2);

showPIC(x0,n,n);
title('reconstructed');

 subplot(3,1,3);
 plot( iterError(1:length(iterError)) );
 
 xcs = x0;
 xbp = xpDyna(:);
 x = samplePIC;
snrValue = SNR(x(:),xcs(:));

saveData('test',xcs,xbp,x,iterError,resultPIC,timeDura,snrValue);


